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Abstract. Bias in image restoration algorithms can hamper further 
analysis, typically when the intensities have a physical meaning of inter¬ 
est, e.g., in medical imaging. We propose to suppress a part of the bias 

- the method bias - while leaving unchanged the other unavoidable part 

- the model bias. Our debiasing technique can be used for any locally 
affine estimator including t\ regularization, anisotropic total-variation 
and some nonlocal filters. 


1 Introduction 

Restoration of an image of interest from its single noisy degraded observation 
necessarily requires imposing some regularity or prior on the solution. Being 
often only crude approximations of the true underlying signal of interest, such 
techniques always introduce a bias towards the prior. However, in general, this is 
not the only source of bias. In many cases, even though the model was perfectly 
accurate, the method would remain biased. This part of the bias often emerges 
from technical reasons, e.g., when approaching an NP-hard problem by an easier 
one (typically, using the t\ convex relaxation of an £q pseudo-norm). 

It is well known that reducing bias is not always favorable in terms of mean 
square error because of the so-called bias-variance trade-off. It is important to 
highlight that a debiasing procedure is expected to re-inject part of the variance, 
therefore increasing the residual noise. Hence, the mean square error is not always 
expected to be improved by such techniques. Debiasing is nevertheless essential 
in applications where the image intensities have a physical sense and critical 
decisions are taken from their values. For instance, the authors of [7] suggest 
using image restoration techniques to estimate a temperature map within a 
tumor tissue for real time automatic surgical intervention. In such applications, 
it is so crucial that the estimated temperature is not biased. A remaining residual 
noise is indeed favorable compared to an uncontrolled bias. 

We introduce a debiasing technique that suppresses the extra bias - the 
method bias - emerging from the choice of the method and leave unchanged the 
bias that is due to the unavoidable choice of the model - the model bias. To 
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that end, we rely on the notion of model subspace essential to carefully define 
different notions of bias. This leads to a mathematical definition of debiasing for 
any locally affine estimators that respect some mild assumptions. 

Interestingly, our debiasing definition for the l\ synthesis (also known as 
LASSO [2D] or Basis Pursuit 0) recovers a well known debiasing scheme called 
refitting that goes back to the “Hybrid LASSO” [D] (see [T5] for more details). 

For the i\ analysis m, including the t\ synthesis but also the anisotropic 
total-variation m, we show that debiasing can be performed with the same 
complexity as the primal-dual algorithm of |4j producing the biased estimate. 

In other cases, e.g., for an affine version of the popular nonlocal-means [2], 
we introduce an iterative scheme that requires only a few run of an algorithm of 
the same complexity as the original one producing the biased estimate. 


2 Background 


We consider observing / = / 0 + w £ R p a corrupted linear observation of 
an unknown signal u o € R w such that /o = <Puo where <T £ R ArxP is a linear 
operator and w is a random vector modeling the noise fluctuations. We assume 
that E[it?] = 0 where E is the expectation operator. The linear operator T> is a 
degrading operator typically with P ^ N and with a non-empty kernel encoding 
some information loss such that the problem becomes ill-posed. 

We focus on estimating the unknown signal uo- Due to the ill-posedness of the 
observation model, we consider variational approaches that attempt to recover 
uq from the single observation / as a solution of the optimization problem 


u*f £ argmin E(u, f) . 

h£ R w 


(1) 


where E : R w xR p —> R is assumed to have at least one minimum. The objective 
E is typically chosen to promote some structure, e.g., smoothness, piece-wise 
constantness, sparsity, etc., that is captured by the so-called model subspace A Pf. 
Providing uj is uniquely defined and differentiable at /, we define C R w as 
the tangent affine subspace at / of the mapping / i —> uj, i.e., 


M} = w}+Im[j;] = {u£R n ; 3z£R p ,u = u) + JJz } 


du* f 

with J * f = W 


( 2 ) 


where JJ is the Jacobian operator at / of the mapping /•->•«} (see [22| for an 
alternative but related definition of model subspace). When uj £ Im[Jj], the 
model subspace restricts to the linear vector subspace _A4j = Im [Jj]. In the rest 
of the paper, u*f is assumed to be differentiable at /o and for almost all /. 


Example 1 The least square estimator constrained to the affine subspace C = 
b + Im[A], b £ R w and A £ ~R Nx Q, is a particular instance of ([!]) where 

E{uJ) = \\<Tu- f\ 2 + i C {n) (3) 


and for any set C, ic is its indicator function: lc{u) = 0 if u £ C, +00 otherwise. 
The solution of minimum Euclidean norm is unique and given by 

u) = b+ A(ffiA) + {f -<Tb) 


(4) 
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where for a matrix M, M + is its Moore-Penrose pseudo-inverse. The affine 
constrained least square restricts the solution u) to the affine model subspace 
M.) = b + ImlA^AY] (as Im[M + ] = lm[M t \). Taking C = R N with for instance 
Q = N, A = Id and b = 0, leads to an unconstrained solution u) =<P + f whose 
model subspace is _A/fj = Im[^*] reducing to 1 N when has full column rank. 

Example 2 The Tikhonov regularization (or Ridge regression) \21\lJfj is an¬ 
other instance of © where, for some parameter A > 0 and matrix T £ M. LxN , 

E(uJ) = ^u-ff + ^\\ru\\ 2 . (5) 

Provided Ker<£nKerP ={0}, u) is uniquely defined as u) = (d> t <P+\r t r)~ 1 d> t f 
which has a linear model subspace given by M) = Im[^*]. 

Example 3 The hard thresholding W' used when <P = Id and fo is supposed to 
be sparse, is a solution of © where, for some parameter A > 0, 

E (n, f) = \\ u ~ /II 2 + y Ho > ( 6 ) 

where ||u||o = ff{i£[P]\ Ui ^ 0} counts the number of non-zero entries of u 
and [P} = {1,... ,P}. The hard thresholding operation writes 

{u})x, = fc f and (m})z= = 0 (7) 

where If = {i£[P] ; | /^ | > A} is the support ofu), I) is the complement of If on 
[P], and for any vector v, vi f is the sub-vector whose elements are indexed bylf. 
As u) is piece-wise differentiable, its model subspace is only defined for almost 
all f as Ai) = {u £ M. N ; u %= 0} =Im[Idz / ], where for any matrix M, Mx f is 
the sub-matrix whose columns are indexed bylf. Note that Idz / £’M. Nx & x f . 

Example 4 The soft thresholding m, used when = Id and fo is supposed to 
be sparse, is another particular solution of © where 

E(u,f) = i||w-/|| 2 + A|M|i , (8) 

with Hr = ]T7 \ut\ the £\ norm of u. The soft thresholding operation writes 

(u})i f = fx f - A sign (f Xf ) and ( u})z= = 0 , (9) 

where If is defined as above, and, as for the hard thresholding: A4) = Im[Idz / ]. 

3 Bias of reconstruction algorithms 

Due to the ill-posedness of our observation model and without any assump¬ 
tions on u o, one cannot ensure the noise variance to be reduced while keeping the 
solution u) unbiased. Recall that the statistical bias is defined as the difference 


Statistical bias = E[it£] — uq . 


(10) 
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An estimator is said unbiased when its statistical bias vanishes. Unfortunately 
the statistical bias is difficult to manipulate when / i —> is non linear. We 
therefore restrict to a definition of bias at /o = <Puo as the error — uq. Note 
that when / i —> u ^ is affine, both definitions match (the expectation being linear). 
Most methods are biased since, without assumptions, uq cannot be guaranteed 
to be in complete accordance with the model subspace, i.e., Uo ^ -M/ 0 - It is 
then important to distinguish techniques that are only biased due to a problem 
of modeling to the ones that are biased due to the method. We then define the 
model bias and the method bias as the quantities 

u* fo -u 0 = u* fo -n M * fo (uo)-n iM * o) ±(uo) , (11) 

Method bias Model bias 

where for any set S, II s denotes the orthogonal projection on S and S - L denotes 
its orthogonal set. We now define a methodically unbiased estimator as follows. 

Definition 1 An estimator Uj is methodically unbiased if 

Vm 0 e r n , u* fo = n M * o (u 0 ) 

We also define the weaker concept of weakly unbiased estimator as follows. 
Definition 2 An estimator Uj is weakly unbiased if 

Vito € A4*f 0 , u*f 0 = ito- 

The quantity u — uo for uq € M*f 0 is called the weak bias of uj at uq- 
Remark that a methodically unbiased estimator is also weakly unbiased. 

Examples. The unconstrained least-square estimator is methodically unbiased 
since u*j q = <P + fo = <£ + ^uo = IIini[^ t ]{ u o) = II M* (uo). Moreover, being linear, 
it becomes statistically unbiased whenever <P has full column rank since <P + <P = 
Id. However the constrained least-square estimator is only weakly unbiased: its 
methodical bias only vanishes when uq G i.e., when there exists to £ 

such that uo = b+A(<I’A) t to. The hard thresholding is also methodically unbiased 
remarking that u*f o is the orthogonal projection on M*f 0 =Im[Idx / ]• Unlike the 
unconstrained least-square estimator, Tikhonov regularization has a non zero 
weak bias. The soft thresholding is also known to be biased El and its weak 
bias is given by —AIdj /o sign(/o)x /o ■ Often, estimators are said to be unbiased 
when they are actually only weakly unbiased. 

4 Definitions of debiasing 

Given an estimate u*j of uo, we define a debiasing of u*j as follows. 

Definition 3 An estimator u*j of uq is a weak debiasina of if it is weakly 
unbiased and A for almost all f, with the model subspace of at 
f. Moreover, it is a methodical debiasina if it is also methodically unbiased. 
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Examples. The unconstrained least square estimator is a methodical debiasing of 
the Tikhonov regularization, since it is a methodically unbiased estimator of uq 
and they share the same model subspace. The hard thresholding is a methodical 
debiasing of the soft thresholding, for the same reasons. 

A good candidate for debiasing u} is the constraint least squares onM^: 

u} = u} + Uf(<I>Uf) + (f — d>u}) € argmin \<Pu — f || 2 (12) 

ueM* f 

where U} gR iVxn with n = rank[Jj] is a matrix whose columns form a basis of 
Im[Jj?]. Let Vj £R nxP be a matrix such that J} = UfVf. The following theorem 
shows that under mild assumptions this choice corresponds to a debiasing of u}. 

Theorem 1 Assume that f i—>• u} is locally affine for almost all f and that <P 
is invertible on A4}. Then u} defined in Eq. m is a weak debiasing ofu}. 


Proof. Since / u} is locally affine, / i-» Uj can be chosen locally constant. 
Deriving ed for almost all / leads to the Jacobian J} of u} given by 


du} du * 


f 

T*\ + 


df du} 


J f = ^r = -df + Wr (ih -■) = J f + 0j}) 


df df 


= U }\7 + U}{$U}) + (ld - $U* f Vf) = 


f 


(13) 


since <TU} has full column rank due to the assumption that ^ is invertible on 
M.}. It follows that 

M) = u) + Im[J;] = u) + U]{$U}) + {f - 3>u}) + lm[U}(<KJ}) + ] (14) 

= u) + Im [U}($U})+] =u}+ Im [U* f ] = M) , (15) 

since <PU} has full column rank. Moreover, for any uq gM.} o , the equation T>u = fo 
has a unique solution u = Uq in M} 0 since <P is invertible on M} 0 ■ Hence, u} o =uq 
is the unique solution of ED. which concludes the proof. □ 


The next proposition shows that the condition a T> invertible on M.}” can be 
dropped when looking at u} and u} through <P. The debiasing becomes further¬ 
more methodical. 


Proposition 1 Assume f u} is locally affine for almost all f. Taking u} 
defined in Eq. ED: then the predictor T>u} of fo = <Tuo is equal to n$M*(f) and 
is a methodical debiasing of d>u}. 

Proof. Since T > U}(T>U}) + = , we have <Pu} = II$m* (/)• As the orthog¬ 

onal projector on its own model space, it is methodically unbiased. Moreover 
Im[JT Im [#£/*]] =Im[$Uf], hence d>u} and <Pu} share the same model subspace. □ 

Remark 1 As an immediate consequence, the debiasing of any locally affine 
denoising algorithm is a methodical debiasing, since d>u} = u}. 
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We focus in the next sections on the debiasing of estimators without explicit 
expression for yVf], meaning that Eq. m cannot be used directly. We first 
introduce an algorithm for the case of t\ analysis relying on the computation of 
the directional derivative Jff- We propose next a general approach, applied to 
an affine nonlocal estimator, that requires J)8 for randomized directions 5. 

5 Debiasing the £ 1 analysis minimization 

From now on, the dependency of all quantities with respect to the observation 
/ will be dropped for the sake of simplicity. Given a linear operator T cR ixAr , 
the t\ analysis minimization reads, for A > 0, as 

E(u,f) = ^ u -f\\ 2 + X\\ru\\ 1 . (16) 

Provided Ker<£ fl Ker W={0}, there exists a solution given implicitly, see [22] , as 

u* = U($U) + f - \U{U t <P t $U)- 1 U t {r t ) x s x (17) 

for almost all / and where 1= {i ; (tit*); ^ 0} C [L\ = {1,..., L} is called the 
co-support of the solution, s = sign(Tit*), U = UJ is a matrix whose columns 
form a basis of Ker[Id X cT] and <TU has full column rank. Note that s x and U are 
locally constant almost everywhere since the co-support is stable with respect 
to small perturbations [22] • It then follows that the model subspace is implicitly 
defined as M* = Im[17] = Ker[IdxcT], and so, the l\ analysis minimization suffers 
from a weak bias equal to —XU(U t <^ t( PU)~ 1 U t (r t ) x s x . Given that u* £ Im[{7] 
and it is locally affine, its weak debiased solution is defined for almost all / as 

u * = U($U) + f . (18) 

The t\ synthesis mu consists in taking r = Id, hence U = ld x , so (ED becomes 

u i = (^i) + / - A((^i) t ^i) _1 si and ufa B = 0 . (19) 

Its model subspace is implicitly defined as M* = Im[Idx], its weak bias is 
—AIdx((^x) t ^x) _1 sx and its weak debiasing is u * = ld x (<P x ) + f. Subsequently, 
taking = Id leads to the soft-thresholding presented earlier. 

The anisotropic Total-Variation (TV) I TS] is a particular instance of (11611 where 
uq £ can be identified to a d-dimensional discrete signal, for which T GR ixW , 
with L = dN, is the concatenation of the discrete gradient operators in each 
canonical directions. In this case I is the set of indexes where the solution has 
discontinuities (non-null gradients) and M* is the space of piece-wise constant 
signals sharing the same discontinuities as the solution. Its weak bias reveals a 
loss of contrast: a shift of intensity on each piece depending on its surrounding 
and the ratio between its perimeter and its area, as shown, e.g., in T9]. Note 
that the so-called staircasing effect of TV regularization is encoded in our frame¬ 
work as a model bias, and is therefore not reduced by our debiasing technique. 
Strategies devoted to the reduction of this effect have been studied in, e.g., [Hy. 
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Since in general u* lias no explicit solutions, it is usually estimated thanks 
to an iterative algorithm that can be expressed as a sequence u k converging to 
u*. The question we address is how to compute u* in practice, i.e., to evaluate 
Eq. (118fl . or more precisely, how to jointly build a sequence u k converging to u*. 

We propose a technique that relies on the observation that, given 1171) . for 
almost all /, the Jacobian J* of u* at / applied to /, leads to Eq. i-e., 

J*[f] = U($U)+f = u* (20) 

since U and si are locally constant [22]. We so define a sequence u k which is, up 
to a slight modification, the closed-form derivation of the primal-dual sequence 
u k of a- Most importantly, we provide a proof of its convergence towards u*. 

Note that other debiasing techniques could be employed for the l\ analysis, 
e.g., using iterative hard-thresholding ra, refitting techniques Ena, post- 
refinement techniques based an Bregman divergences and nonlinear inverse scale 
spaces I17I3I24I or with ideal spectral filtering in the analysis sense [| I71 . 


5.1 Primal-dual algorithm 

Before stating our main result, let us recall some of the properties of primal- 
dual techniques. Dualizing the l\ analysis norm u A||Tn||i, the primal problem 
can be reformulated as the following saddle-point problem 

z* = argrnax min -\\$u - f\\ 2 + (Aj, z) - i Bx {z) (21) 

zeR i u 6R W 2 

where z* £ K L is the dual variable, and B\ = {z ; ||^||oo ^ A} is the too ball. 


First order primal-dual optimization. Taking err < , 9 £ [0,1] and initializ¬ 

ing (for instance,) u° = v° = 0 £ M. N , z° = 0 £ K L , the primal-dual algorithm 
of [1] applied to problem (12T1) reads 


z k+1 = n Bx (z k + crr v k ), 

u k + 1 = (Id + T^)" 1 (u k + r(4> 4 / - r\z k+1 ))) , 
v k +l _ yk+ 1 _|_ Q(yk+1 - u fe) ) 


where the projection of z over B\ is done component-wise as 


n Bx (z)i 


Zi if \zi\ ^ A, 

A sign(zj) otherwise. 


( 22 ) 


(23) 


The primal-dual sequence u k converges to a solution u * of m We assumed 
here that u* verifies m with <FU full-column rank. This could be enforced as 
shown in [22], but it did not seem to be necessary in our experiments. 


5.2 Debiasing algorithm 

As pointed out earlier, the debiasing of u* consists in applying the Jacobian 
matrix J* at / to / itself. This idea leads to the proposed debiasing algorithm 
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that constructs a sequence of debiased iterates from the original biased primal- 
dual sequence with initialization u° = v° = 0 € R N , z° = 0 € R L as follows 


where 


~z k+1 = n zk+arvk {z k + an k ), 

h k+1 = (Id + T ^)- 1 ( u k + T(&f - r 4 i fe+1 )) , 

v k+1 = u k+1 + 9(u k+1 - u k ), 


Hz k +aTv k {Zi) 


Zi if \z k + arv\ < A + /3, 
0 otherwise. 


(24) 


with f3 ^ 0. Note that when /3 = 0, deriving z k , u k and v k for almost all / at 
/ in the direction / using the chain rule leads to the sequences z k , u k and v k 
respectively (see also jjj]). However, as shown in Theorem EJ it is important to 
choose /? > 0 to guarantee the convergence of the sequencqj- 

Theorem 2 Let a > 0 be the minimum non zero wjZk<H of for all i £ [L\. 

Choose /? such that aa > 0 > 0. The sequence u k defined in (El) converges to 
the debiasing u * of u*. 


Before turning to the proof of this theorem, let us introduce a first lemma. 

Lemma 1 The debiasing u* of u * is the solution of the saddle-point problem 

min max \\<Tu - f\\ 2 + ( Tu , z) - t F Jz), (25) 

ueR N JeR 1 - 

where l Fx is the indicator function of the convex set F x = {p £ ; p x = 0} . 

Proof. As <PU has full column rank, the debiased solution is the unique solution 
of the constrained least square estimation problem 

u * = U{T>U) + f = argmin \\<Pu — f || 2 . (26) 

u eu* 

Remark that u G U* = Ker[IdjoT] (Tu) x= = 0 L Fxc (ru) = 0, where 

F IC = {p e ; p X c= 0}. 

Using Fenchel transform, i Fxc {ru) = max: (ru,z) — i Fxc (z), where i F is 
the convex conjugate of i F a . Observing that i Fx =i* Fxc concludes the proof. □ 

Given LemmaU replacing n zk+aFvk in (El) by the projection onto Fx, i.e., 

II Fx (z) X c = z X c and n Fx (z) x = 0, (27) 

leads to the primal-dual algorithm of [¥j applied to problem (El) which converges 
to the debiased estimator u*. It remains to prove that the projection LI z k +aFvk 
defined in (El) converges to 77 Fx in finite time. 

Proof (Theorem @). First consider i 6l, i.e., |7 n rt*| i >0. By assumption on a , 
\ru*\i>a > 0. Necessary z* = Asign(7 n M*) i in order to maximize El) . Hence, 

| z* + oFu*\ i > A + era. Using the triangle inequality shows that 

A + aa < \z* + aru*\ t < \z* - z k \ l + a\r U * - rv\ + \z k + aFv k \ i . (28) 

1 In practice, (5 can be chosen as the smallest positive floating number. 

2 If IFu*^ =0 for all i £ [L\, the result remains true for any a > 0. 
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Choose e > 0 sufficiently small such that era — e(l + a) > j3. From the conver¬ 
gence of the primal-dual algorithm of [3], the sequence (z k ,u k ,v k ) converges to 
(z*, u*,u*). Therefore, for k large enough, | z* — z k \ i <£, \Tu* — Fv k |, <e, and 

\z k +arv k \i ^ X + aa-e(l+a) > X + /3 . (29) 

Next consider i£l c , i.e., |r , u*|i = 0, where by definition \z*\i ^ A. Using again 
the triangle inequality shows that 

\z k + °rv k I, < |z fc - A + a\rv k - ru*u + |z*|,. (so) 

Choose £ > 0 sufficiently small such that e(l + a) < /3. As ( z k ,u k ,v k ) —> 
(z*,u*,u*), for k large enough, \z k - z* |j<£, | rv k — ru*\ i <£ , and 

\z k +arv k \ i <X + £(l+a) ^ X + /3 . (31) 

It follows that for k sufficiently large | z k + cr T v k | { ^ A + /3 if and only if i £l c , and 
hence II z k +aKv k(z) = IIf x {z). As a result, all subsequent iterations of will 
solve and hence from Lemma [1] this concludes the proof of the theorem. □ 

6 Debiasing other affine estimators 

In most cases, U* cannot be computed in reasonable memory load and/or 
time, such that Eq. (ED cannot be used directly. However, the directional deriva¬ 
tive, i.e., the application of J* to a direction 8, can in general be obtained with 
an algorithm of the same complexity as the one providing u* . If one can compute 
the directional derivatives for any direction, a general iterative algorithm for the 
computation of u * can be derived as given in Algorithm 1. 

The proposed technique relies on the fact that given n = dim(Al*) uniformly 
random directions <5i,..., S n on the unit sphere of R p , J*<5i,..., J*S n forms a 
basis of Im[J*] almost surely. Given this basis, the debiased solution can so be 
retrieved from ED- Unfortunately, computing the image of the usually large 
number n of random directions can be computationally prohibitive. 

The idea is to approach the debiased solution by retrieving only a low dimen¬ 
sional subspace of A4* leading to a small approximation error. Our greedy heuris¬ 
tic is to chose random perturbations around the current residual (the strength of 
the perturbation being controlled by a parameter e). As soon as e > 0, the algo¬ 
rithm converges in n iterations as explained above. But, by focusing in directions 
guided by the current residual, the algorithm refines in priority the directions 
for which the current debiasing gets significantly away from the data /, i.e., di¬ 
rections that encodes potential remaining bias. Hence, the debiasing can be very 
effective even though a small number of such directions has been explored. We 
notice in our experiments that with a small value of e, this strategy leads indeed 
to a satisfying debiasing, close to convergence, reached in a few iterations. 

The nonlocal-means example. The block-wise nonlocal-means proposed in [5] 
can be rewritten as an instance of the minimization problem © with 

E (U, f)=\j2 W iA V i U ~ V lff with w iJ = V ( *V J/I ) < ' 32 ^ 

i,j ^ ' 
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Algorithm 1 General debiasing pseudo-algorithm for the computation of u*. 


Inputs: / £ R p , u* £ R w , 5 £ R p -¥ J*S £ R N , e > 0. 

Outputs: u * € R^ and U £ R iVX71 with n'^naa orthonormal family of Im[J*] 


Initialize U £ 
repeat until u 
Generate 
Compute 
Compute 
Update 
Update 


reaches convergence 

5 <— ?7/||?7||, p ~ A/p(0, Id) (perturbation ensuring convergence) 

u <— J* (/ — <fru k + eS) (perturbed image of the current residual) 
e <— u — C/(t/*w / ) (projection on the orthogonal of the current At*) 
U<-[U e/||e||] 

u* 4 — u* + U({<PU)+(f - <2hV)) 


where i £ [ni] x [ 712 ] spans the whole image domain, j — i £ [—s, s] x [—s, s] spans 
a limited search window domain and a 2 is the noise variance. We denote by Vj 
the linear operator extracting the patch at pixel i of size (2 p + 1) x (2 p + 1). 
Note that we assume periodical conditions such that all quantities remain inside 
the image domain. The kernel ip : R + —>• [0,1] is a decreasing function which is 
typically a decay exponential function. Taking ip piece-wise constant^, leads to 
computing u* and its Jacobian at / applied to 6 for almost all / as follow 

* , /• T *c\ Wj^Sj _ V-^ . . 

Uj = _ - and (J S)i = _ with w itj = } Wi- k ,j-k (33) 

where k £ [—p,p] x [—p,p] spans the patch domain. Note that the values of w 
and w can be obtained by discrete convolutions leading to an algorithm with 
complexity in 0(Ns 2 ), independent of the half patch size p. 

With such a choice of tp, the block-wise nonlocal filter becomes a piece-wise 
affine mapping of / and hence Algorithm 1 applies. 

7 Numerical experiments and results 

Figure |T] gives an illustration of TV used for denoising a ID piece-wise con¬ 
stant signal in [0,192] and damaged by additive white Gaussian noise (AWGN) 
with a standard deviation a = 10. Even though TV has perfectly retrieved the 
support of V«o with one more extra jump, the intensities of some regions are 
biased. Our debiasing is as expected unbiased for every region. 

Figure [2] gives an illustration of our debiasing of 2D anisotropic TV used for 
the restoration of an 8 bits approximately piece-wise constant image damaged by 
AWGN with er = 20. The observation operator is a Gaussian convolution kernel 
of bandwidth 2px. TV introduced a significant loss of contrast, typically for the 
thin contours of the drawing, which are re-enhanced by our debiased result. 

Figure [3] gives an illustration of our iterative debiasing for the block-wise 
nonlocal-means algorithm used in a denoising problem for an 8 bits image enjoy¬ 
ing many repetitive patterns and damaged by AWGN with a = 20. Convergence 
has been considered as reached after 4 iterations only. Our debiasing provides 
favorable results with many enhanced details compared to the biased result. 

For instance, by quantification on a subset of predefined values in [0,1]. 
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Fig. 1. Solutions of ID-TV and our debiasing on a piece-wise constant signal. 



PSNR 19.13 / SSIM 0.76 PSNR 20.61 / SSIM 0.80 PSNR 21.90 / SSIM 0.87 


Fig. 2. (left) Blurry image f = 4>uo+w, (center) TV u *, (right) debiased u*. 



PSNR 22.14 / SSIM 0.52 PSNR 27.89 / SSIM 0.82 PSNR 29.17 / SSIM 0.87 


Fig. 3. (left) Noisy image f = uo+w, (center) nonlocal-means u*, (right) debiased u*. 

8 Conclusion 

We have introduced in this paper a mathematical definition of debiasing 
which has led to an effective debiasing technique that can remove the method 
bias that does not arise from the unavoidable choice of the model. This debiasing 
technique simply consists in applying a least-square estimation constrained to 
the model subspace chosen implicitly by the original biased algorithm. Numerical 
experiments have demonstrated the efficiency of our technique in retrieving the 
correct intensities while respecting the structure of the original model subspace. 
Our technique is nevertheless limited to locally affine estimators. Isotropic total 
variation, structured sparsity or nonlocal-means with smooth kernels are not yet 
handled by our debiasing technique, and left for future work. 
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